ECT*-04-07 



Forward-backward equations for nonlinear propagation in axially- invariant optical 

systems 

Albert Ferrando and Mario Zacares. 
Departament d'Optica, Universitat de Valencia. Dr. Moliner, 50. E-^6100 Burjassot (Valencia), Spain. 



Pedro Fernandez de Cordoba. 

Departamento de Matemdtica Aplicada, Universidad Politecnica de Valencia. Camino de Vera, s/n. E-^6022 Valencia, 



ECT* 



Mediterranean 



Daniele Binosi. 

Villa Tambosi, Strada delle Tabarelle 286, 1-38050 Villazzano (Trento), Italy. 



Alvaro Montero. 

of Science and Technology. Camino de Vera, s/n. E-46022 Valencia, 
(Dated: 2nd February 2008) 



We present a novel general framework to deal with forward and backward components of the 
electromagnetic field in axially-invariant nonlinear optical systems, which include those having any 
type of linear or nonlinear transverse inhomogeneities. With a minimum amount of approximations, 
we obtain a system of two first-order equations for forward and backward components explicitly 
showing the nonlinear couplings among them. The modal approach used allows for an effective 
reduction of the dimensionality of the original problem from 3 + 1 (three spatial dimensions plus one 
time dimension) to 1 + 1 (one spatial dimension plus one frequency dimension). The new equations 
can be written in a spinor Dirac-like form, out of which conserved quantities can be calculated in 
an elegant manner. Finally, these new equations inherently incorporate spatio-temporal couplings, 
so that they can be easily particularized to deal with purely temporal or purely spatial effects. 
Nonlinear forward pulse propagation and non-paraxial evolution of spatial structures are analyzed 
as examples. 

PACS numbers: 42.65.-k, 42.65. Sf, 52.35.Mw 



INTRODUCTION 



Nonlinear propagation of light pulses in dielectric me- 
dia such as optical fibers has been traditionally mod- 
eled using the Nonlinear Schrodinger Equation (NLSE) 
Q. However, it is well known that NLSE needs mod- 
ifications to describe a number of higher-order nonlin- 
ear effects which become important at increasing powers 
and for short pulses. Recently, the access to new optical 
systems in which nonlinearities can be considerably en- 
hanced together with the experimental availability of ul- 
trashort pulses push the description based on the NLSE 
and its modified versions to a limit. A typical exam- 
ple of this new scenario is provided by the phenomenon 
of supercontinuum generation in photonic crystal fibers 
0, which requires a specific modeling that goes beyond 
approaches based on conventional versions of the NLSE 
These new approaches are expressed in new evolu- 
tion equations that differ from the NLSE in the amount 
of approximations needed to achieve them. We can men- 
tion the so-called generalized NLSE 0, 0] , the nonlinear 
envelope equation (NEE) Q, the forward Maxwell equa- 
tion (FME) and the unidirectional pulse propagation 
equation (UPPE) Q . Briefly, the aim of these equations 
is to describe pulse propagation in the regime where the 
frequency width of the pulse is comparable to the carrier 



frequency, which, in turn, translates into the fact that 
usual approaches suchas the slowly varying approxima- 
tion no longer hold. The specific form of these equations 
is, on the one hand the result of applying some other dif- 
ferent approximations, e.g, assuming propagation in an 
homogeneous medium (NEE and UPPE) or single-mode 
propagation in fibers (generalized NLSE and FME). On 
the other hand, two common features of all of them are 
the neglect of the backward components of the electro- 
magnetic field and their first-order character. 

The role played by backward components has been pre- 
viously analized for an homogeneous medium Q. In this 
paper we explicitly unveil their role but in the more gen- 
eral case of an axially-invariant inhomogenous nonlinear 
medium by explicitly finding the coupled first-order equa- 
tions that drive the forward and backward components 
of the electromagnetic field in an axially-invariant nonlin- 
ear system. We will show that these first order forward- 
backward equations (FBEs) are equivalent to the orig- 
inal second-order equations for the electric components 
of the electromagnetic field. The FBEs will provide us 
with a common framework that can encompass different 
nonlinear phenomena. In fact, since these equations ex- 
plicitly manifest the couplings between spatial and time- 
frequency degrees of freedom typical of spatial-temporal 
phenomena, they can be easily particularized to describe 
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either purely temporal or purely spatial effects within the 
same framework, revealing the total generality of this for- 
malism. 

The paper is organized as follows. In section |H] we 
derive the most general modal second-order equation for 
nonlinear propagation in an axially-invariant inhomoge- 
nous medium and we explain the nature of the only ap- 
proximation needed to obtain it. In section lllll we demon- 
strate the equivalence between the modal second-order 
equation and the first-order FBEs. In section|^we intro- 
duce a spinor representation to obtain a Dirac-like form 
of the FBEs. In section^we derive the conserved quanti- 
ties associated to the FBEs and analize them in the light 
of phase symmetries. Finally, in section we examine 
two different nonlinear phenomena occurring in axially- 
invariant inhomogeneous dielectric media in the light of 
the FBEs: forward pulse propagation and non-paraxial 
evolution of spatial structures. 

II. MODAL SECOND-ORDER EQUATION 

The most general equation for the propagation of 
the electric components of a electromagnetic field in a 
inhomogeneous, isotropic and spatially local dielectric 
medium is given by 

V 2 E W - V(V • E w ) + fc^E^ = -G^ L (E), (1) 

where E w = E w (x,?/,z) is the complex u frequency 
component of the real electric field E(ar, y, z, t) — 
E w (x, y, z) exp(— icut), V is the three-dimensional spa- 
tial gradient operator, ko is the vacuum wavenumber, 
no = no(x,y, z) is the refractive index profile of the di- 
electric medium and G^ L is the function that represents 
the local nonlinear response of the medium to the propa- 
gating field. In the most general case, the relative dielec- 
tric constant (e = Uq) and the nonlinear function G^ L 
can be either spatially local (the case we are considering) 
or non-local if there exist some type of spatially-delayed 
response effects. Besides the general wave equation 
Maxwell's equations require the constraint V • D w = to 
be satisfied. Through the constitutive relations, the dis- 
placement field D w is itself a function of the propagating 
electric field E w . Assuming a spatially local response, 
D w = e(E w )E aJ . In the case the system presented some 
type of anisotropy, e(E w ) would be a second-order tensor. 
Here we consider an isotropic medium, so that e(E w ) will 
be a scalar function, although the generalization to an 
anisotropic medium is straightforward. The constraint 
relation has thus the form: 

V • (e(E a ,)E w ) = 0. (2) 

In most cases in nonlinear optics, the general constraint 
<J2| is replaced by the approximated, and much simpler, 
"scalar" condition V • (e(EjEj « e(E w )V • E w = 0, 
implying that V ■ E u k 0. This approximation is known 
to work well in a panoply of nonlinear effects with few and 



remarkable exceptions, related mainly to extreme self- 
focussing eventsjl]. In order to simplify our approach, 
this approximation will be assumed throughout, although 
our main results will remain valid in a slightlier general 
context in which the condition V • E w « will not be 
strictly necessary. 

Our interest lies on describing propagation in an 
axially-invariant system, that is, that in which both the 
refractive index profile no(x,y,z) = no(x,y) as well as 
all other macroscopic nonlinear structural functions, such 
as nonlinear susceptibilities, are z-independent. In other 
words, we will focus on systems for which the axial z- 
dependence is carried by the propagating field ~E(x, y, z) 
exclusively. Mathematically, this property is reflected 
into the fact that there is no explicit dependence on z in 
the nolinear function G NL . In these systems, the wave 
equation and the constraint V • E w w merge into a 
single equation: 

+ V? + klnlix, y)j E w = -G^ L (E). (3) 

Where now V* stands for the transverse gradient opera- 
tor and the axial second-order derivative has been made 
explicit. Let us now perform a modal expansion of the 
electric field in terms of the eigenmodes of the linear z- 
independent system, described by the linear transverse 
operator L = Vf + k^n^x, y). These modes fulfill the 
linear eigenvalue equation: 

L^ l {x. : y)=(5l{uj)^ n {x 1 y). (4) 

Since Lq does not mix the spatial (polarization) compo- 
nents of the electric field, it is proportional to the iden- 
tity operator in polarization space. For this reason, its 
eigenmodes have a three-fold degeneracy in polarization 
indices. Therefore, every multiplet of eigenmodes is con- 
stituted by three linearly independent modes, each one 
proportional to a three-dimensional vector belonging to 
a basis of M 3 . For simplicity, we can consider this ba- 
sis to be the canonical one ^^,Jx,y) = <fr"(:r, y)u( CT ) 
(a = 1, 2, 3) where the components of the canonical basis 
{u(i),U( 2 ),U(3)} satisfy u (t7)a = S aa (a = 1,2,3). If the 
system is lossless then no = Uq, and Lq is self-adjoint; if 
it is not, we consider the real part of no in Eq.l@J so that 
Lq is again self-adjoint. In this way, the set of Lq eigen- 
modes f° rm an orthogonal basis of functions of 
the transverse coordinates that can be used to expand 
the electric field Fi u (x, y, z) at a given z: 

E u (x,y,z) = ^c„ !CT (w;z)^(x,y)u ((T) , 

n,<7 

or, in components, 

E^{x,y,z) = Y J Cna(w,z)^ n (x,y). (5) 

n 

Of course, the modal expansion © has to be under- 
stood as a generalized form of expansion over the entire 
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spectrum of the Lq operator, in which both discrete and 
continuum parts of the spectrum can coexist. In prac- 
tice however, when dealing with numerical applications, 
the continuum part is discretized by means of some con- 
venient election of boundary conditions, so that the dis- 
cretized form of the expansion (jjj is, in fact, the one 
that is used. Since the election of boundary conditions is 
an artifact of the simulation, one has to make sure that 
physical results are independent of them. 

Now, we substitute the modal expansion © into 
Eq.©, we multiply this equation by & m * , taking into 
account that LqQ^ — /3„$^, and after integrating the re- 
sult over the entire transverse space (considering, without 
loss of generality, that the functions are orthonormal: 
J R2 $m < i ) 5i = we obtain the evolution equation for 

the modal expansion coefficients c: 

+ c na (io- z) = -F^(c). (6) 

The equation above represents an effective lower- 
dimensional equation (1 + 1 dimensions, instead of the 
3 + 1 dimensions of the original equation ©) for the 
coefficients c nr7 (uj;z). Physically, they represent the fre- 
quency, modal and polarization spectral content of the 
propagating electric field. The nonlinear function F£ a {c) 
is given by: 

F^(c)=[ CC(E^)- (7) 

This function is an input of the effective equation © 
because it is known once the linear amplitudes are de- 
termined out of the linear eigenvalue problem in Eq. J3J . 
Same applies to the propagation constants (3 n {u)), which 
also appear as an input. 

As an example of nonlinear function Ffi a (c), it is in- 
structive to consider the case of a Kerr nonlinearity, 
for which = u 2 /c 2 xS/^ T (w, wa, w 3 , lu + u 2 - 
uj 3 )E^f*E^fE^ + ^- UJs . From the definition of the non- 
linear function 0, one finds that F has the form: 

KM = E E X 

m'n'm r'a'r 

Ju>2^3 a r' a' r 

Thus, for a Kerr nonlinearity, all the information about 
the nonlinear properties of the system is encoded in the 
tensor function T. As a first approximation, one could 
say that this tensor function depends on frequency both 
through the susceptibility x^ 3 ' and through the overlap- 
ping integral of linear mode amplitudes: 

T nm i „i m (<JJ, W3j UJ + UJ2 — UJ3) ~ UJ 2 /c 2 
a r' a' r 

Besides, the dependence on modal indices of the T ten- 
sor is a result of the overlapping of spatial amplitudes, 



whereas polarization indices are provided by the third- 
order susceptibility. 

We do not want to get into much detail on the par- 
ticular form of the nonlinear function since, for our pur- 
poses, it is not necessary to provide an explicit form for 
F in Eq.©. Only general properties of F will be used 
and the latter can be inferred from the construction pre- 
viously described. Along the same line of reasoning, it 
should be added that the effective Eq.© remains valid in 
some cases when the "scalar" approximation V • E w w 
no longer holds. Less restrictive approximations can be 
made in which V ■ D w « but V • E w ^ and, neverthe- 
less, the effective Eq.© retains, formally, its validity. It 
would differ in the fact that the values of (3 would be now 
calculated including the missing vector term in the linear 
eigenvalue equation QJ and that the nonlinear function 
F would include also terms of vectorial origin. Again, the 
analysis of such topics is out of the scope of this paper, 
for our aim is to transform the almost completely general 
second-order equation © into a first-order formalism in 
which forward and backward components explicitly ap- 
pear. 



III. DERIVATION OF FORWARD-BACKWARD 
EQUATIONS 

A convenient way of considering the effective equa- 
tion © is as a limit of a equation for a frequency- 
discretized spectral function c 3 n(T (z) = c nr7 (ujj;z). This 
is a typical situation in numerical simulations in which 
frequency appears discretized in a Fourier series and one 
approaches the continuum-frequency limit numerically. 
Or, likewise in some experimental cases where one phys- 
ically works with frequency combs instead of continuum 
sources. In such situations, the frequency discretized ver- 
sion of Eq.© is: 

(■^ + Pni) C 'nAz) = -FUc)- (8) 

The general second-order equation ©, displaying all 
frequency, modal and polarization indices, will be our 
starting point. The function F(c) includes all nonlinear 
contributions in the frequency-domain constructed out 
of Maxwell's equations according to the procedure de- 
scribed in the previous section. As seen above, for a Kerr 
nonlinearity in Maxwell's equation, this function would 
provide a cubic behavior in the spectral function c, which 
in the discretized frequency version would be: 

E E E>< 

m'n'm r'a'r k'j' 

k'* ( rpjk'l(j+k' I fj+k'-j' (Q\ 
^m'r' I „ m'n'm J IlVSlT ' J 

\ cr r' a' r / 

where T is a tensor in frequency, modal and polariza- 
tion indices that can be explicitly calculated out of the 
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amplitudes of the linear fiber modes. However, in or- 
der to obtain the first order counterpart of the general 
second-order equation JHJ, it is not necessary to specify 
a particular form for the nonlinearity function F. It is 
interesting to stress again that both the modal disper- 
sion relations (3 n (u>) (or its discretized-frequency version 
flnj) and F are inputs that are known once the mode 
linear problem is solved. In order to simplify the nota- 
tion, we will incorporate the polarization index into the 
modal one when dealing with the spectral function c and 
the nonlinear function F. This is possible because polar- 
ization and modal indices always come together in pairs. 
So that, from now on and without loss of generality, n 
represents the index pair (n, a) . 

We perform now an axial Fourier transform on the 
spectral functions c in Eq.ljHJ, defined as 



i 

2^ 



or, inversely, 



dzc' n (z)e~ ll3z . 



With this definition, (3 — > —id/dz (or, d/dz — » i/3), 
and, thus, the axial Fourier transfom of Eq.© is given 
by 



(10) 



where F(c) is the axial Fourier transform of the nonlinear 
term in Eq.©: F = Fp[F]. Dividing this equation by 
— 1 + (we assume some kind of "ie" prescription to 
deal with the poles of the inverse of this function in the 
P plane), we obtain 



1 



P - ft 



(ii) 



The function preceding the nonlinear factor in the 
above equation is the axial Green function in /3-space 
and has single poles at +/3 n j and —@ n j, corresponding to 
forward and backward propagation, respectively, for the 
positive frequency part of the spectral function c (peaked 
at a positive frequency +u>q). Clearly, we can decompose 
the Green function in the contributions corresponding to 
the two poles by means of the following identity: 



'■(■ 1 



1 



which allows us to write En. lfTT)) as 

c = c F + cb, 



(12) 



where 



(c F y r 



i i 



1 



1 



2/3 nj + Pnj 



or, equivalently, 

08 - nj ) (c F ) J n = 

w+pn 3 )(c B y n = 



Now, we take the inverse axial Fourier transform of 
the above equations, taking into account the previous 
definitions. So that, considering that ft — > —id/dz, 

T^icF.B) = T~ 1 (J 7 p(cf,b)) = cf,b, and F' 1 ^) = 
J-~ l {J- 13(F)) = F, we can write the following two first- 
order equations for forward (F) and backward (B) com- 
ponents: 



2/3, 



-Fl(c), 



(13) 



(14) 



where the total spectrum is given by the sum of forward 
and backward contributions 



cf + cb ■ 



(15) 



Eqs. l(T3ll and ifl^Ijl are the exact first-order equations 
equivalent of the general second-order equation |JHJ . 

Restoring the continuum frequency notation (e 7 — > 
c(u>)), we obtain the first-order forward-backward equa- 
tions in the continuum-frequency limit: 



d 

-i-^-Pn{u) ) c F (u,z) 



2/3„(w) 



whereas the total spectrum is 



2/3„(w 



F n (c(u,z)) (16) 



■F n (c(w,z)) (17) 



c(lo, z) = c f {lo, z) + Cb(w, z). 

Notice that if the backward spectrum is sufficiently 
small and it can be approximately neglected (eg ~ 0), 
then c ~ cf and only Eqs. l(T3jl or (tlfill matter. It is also 
interesting to notice that small values of the backward 
spectrum imply small nonlinearities, as it can be checked 
by going to the cb — > limit in Eas.ltl4ll or <fl7t . In 
such a limit, nonlinearities dissapear independently of 
the nature of their origin: F(c) — > 0. We will return to 
this point in the last section of this paper. 
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IV. SPINOR REPRESENTATION OF THE 
FORWARD-BACKWARD EQUATIONS 

It is possible to put the first order FBEs in a more 
compact form, by using a spinor representation for the 
forward and backward spectral functions, cf and cb- We 
start noticing that the FBEs (with discrete-frequency in- 
dices) can be rewritten as, 



(- 



dz 



< d 



n'j' 



(18) 

We have assumed that the nonlinear function F can 
be expressed as the action of a field-dependent opera- 
tor M(c) on the spectral function c; that is, F^(c) = 

Yln'j' Mnn' (c)^' • This is the situation that applies 
to all type of nonlinearities that can be expanded in 
a power series. The simplest case is the aforemen- 
tioned Kerr cubic nonlinearity, for which M^ n , (c) = 

EL^^S'f'" 3 '''^^ 1 '' Moreover, in conser- 
vative systems, that is, those for which the Hamiltonian 
is conserved and real, the M operator is forced to be self- 
adjoint. The N operator appearing in Eo. lfT8|l is nothing 

but JV£,( C ) ee 1/(2^-)m£,(c). ' 

The equations above can also be written in a matrix 
form, 



—i 




Then, introducing the bi-spinor if), defined as 



3 = 1 F 



(19) 



(20) 



and the Pauli matrices, 



0"! 



1 



1 \ ( -i 

, °2 = 



i 



0-3 = 



1 

-1 



noticing that 



1 1 

-1 -1 



= ct 3 + ia 2 , 



we see that Eq. iTLQ)) can be written in a bi-spinor form 
notation (a sum over repeated indices is assumed) 



(21) 



Certainly, TV is a non-trivial operator in frequency and 
mode spaces. However, it is proportional to the iden- 
tity operator when acting on the F — B internal degrees 
of freedom of the bi-spinor ip. If we use continuum- 
frequency notation, Eci J21Jl can be written as 

d 

^ J du'N nn >(u>,w';if))((T3 + i<T 2 )^n<(u ,z). (22) 

n' 

Despite its different appearence, the first-order spinor 
equation for the pulse spectrum contains exactly the 
same information on dynamics than the original second- 
order equation @- The spinor representation of the 
FBEs has common features with the Dirac equation for 
a particle in 1 + 1 dimensions. In 1 + 1 dimensions, the 
original algebra of 16 4 x 4 matrices (Dirac matrices) of 
the 4D Dirac equation is reduced to 4 2 x 2 matrices con- 
stituted by the identity and Pauli matrices. The FBEs 
spinor l(20j) , formed by the forward and backward spectral 
functions, plays the role of the Dirac spinor, constituted 
by the positive and negative energy components of the 
particle wave-function. 



V. CONSERVED QUANTITIES 

The goal of this section is to find the conserved quan- 
tity associated to FBEs in the most general case. Our 
only assumption will be the conservative and lossless 
character of the system, which mathematically will be 
reflected in the fact that the M operator is self-adjoint 
(M = M') and that the dispersion relation function j3 
is real for all modes. We consider as a starting point 
the discrete-frequency FBEs in their spinor representa- 
tion l(2ll . Now we proceed to redefine the labeling of this 
equation following the same procedure we used to incor- 
porate polarization indices into modal indices in Section 
nil Since frequency (J) and modal indices (n) are always 
paired together in Eq.lj2lJ, we can incorporate both in a 
new single index: (n,j) — * p, (n',j') — > p\ In this way, 
Eq. lPU becomes (a sum over repeated indices is assumed, 
as before): 



- = [{Pp'&PP' + n pp') a 3 + iN ppl a 2 ] Vy • (23) 

This form of the equation enables us to introduce the 
matrix operators ^, B and N, whose elements are given 
by ip p , /3p'Sp P ' and N pp r, respectively. Therefore, we can 
write (\i> ee dW/dz): 

- = | ^B + ^B^M^j cr 3 + ^B^Mo-a j (24) 

where we have made use of the relation between 
the N and M operators: N pp > = l/(2(3 p )M pp/ 
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N = (1/2)B- X M. 
above equation: 



We will also need the adjoint of the 



B + -MB" 



03 - j-MB" 1 ^ 



(25) 



where we have used the properties that B = B''' (since B is 
diagonal and (3 is real), M = M 1 *, together with the self- 
adjointness of the Pauli matrices. Notice that both B and 
M are proportional to the identity matrix when they act 
on the F — B components of the bi-spinor *Sf, and for this 
reason they commute with Pauli matrices. Next, we left- 
multiply Eq.JSJ by ^Ba 3 and right-multiply Eq. jSH 
by Ba 3 ^, to obtain (erf = 1, 02(73 = *0"i): 



i¥Ba 3 y = * f \ B 



1 



-M- 



-i^Baa* = * f \ B 2 + -M 



1 



-Mcri }> * 



We achieve the desired result by substracting the pre- 
vious equations: 

i ^Ba 3 ^! + ^Ba 3 ^ = ^ ^- (^Ba 3 ^) = 0. 

(26) 

The conserved quantity is thus Q = ^Ba 3 ^f, or, after 
reintroducing indices, Q = Yl n ,j Ynj^nj^Z^nj (discrete 
frequency) or Q = J2 n J d^l(uj)P n (uj)a 3 ip n (u;) (contin- 
uous frequency). In terms of the original forward and 
backward components of the bi-spinor ip, the conserved 
quantity has the following form: 



(27) 



/3 n (uj) (c* Fn (uj, z)c Fn (uj, z) - c* Bn (uj, z)c Bn {w, z)) . 

The physical meaning of this conserved quantity can 
help us to understand its particular form. Eci, l|27H is the 
modal frequency version of the axial component of the 
Poynting vector, as one could check by reminding that, 
for the case of Maxwell's equations in the scalar approxi- 
mation, V z ~ / K 2 i4>*d z 4>- This represents the amount of 
electromagnetic energy traversing a section of the system 
per unit time. For this reason, since we are dealing in fact 
with the axial flux of the electromagnetic field, it is nat- 
ural that the quantity Q, the total electromagnetic axial 
flux, can be considered as the sum of the positive for- 
ward axial flux (+Qf) and the negative backward axial 
flux (-Qb)- Q = Qf-Qb, where Q F = J2n I PnC* Fn c Fn 
and Q B = J2 n I PnC* Bn c Bn . 



VI. PARTICULAR CASES 

An interesting case to consider is that occuring when 
one neglects all F — B interactions. In general, the FBEs 



llKit and II17|I mix F — B components due to the presence 
of the nonlinear function F(c) — F(c F + c B ) in the right 
hand side of both equations. When F — B interactions 
can be neglected, what happens when either c F or cb 
are very small, FBEs decouple and we obtain two sep- 
arate equations for c F and c B . In such a case, forward 
and backward axial fluxes are conserved independently: 
dQ F /dz — and dQ B /dz = 0. The demonstration of 
this property closely follows the general proof previously 
described for the total flux Q. Let us consider FBEs in 
their discrete-frequency form (EgS- flUt and (Q3J) when 
all F — B mixing terms are neglected; that is, when 
F(c) « F(c F ) in Ea. lll3Jl because c ~ c F (cb ~ 0) or 
F(c) k, F(cb) in Eq.O because c « c B (c F w 0). We 
will consider forward components only (backward analy- 
sis is completely analogous) . By relabelling the frequency 
and modal indices together into a new index and by in- 
troducing the matrix notation in the same way we used 
before to obtain Eq.{23J, the equation for forward com- 
ponents adopts the form: 



ZCf 



B + ^B-^cp)) c P . 



(28) 



In this case, the equation for the backward component 
would correspond to a very weak field (c B ~ 0) and, 
thus, it would be basically linear: — zcb ~ — Bcb. For 
our purposes, we also need the adjoint equation of J2HJ: 



1 c 



t - „„t 



1 



c F f B + -M(c F )B" 1 



Left-multiplying by cp^B and right-multiplying 

the previous equation by Bcf, one gets: 



zcf^Bcf = cf^ 



B 



-M(c F ) c F 



ictBcF = Cf^ ( B 



2 M ( c f) I <>■• 



The desired conservation law comes now after substract- 
ing the above equations: 



^(c F Bc F ) = CF f BCF 



clBcp = 0. 



Certainly, Q F = c f Bcf = J2 n I PnC Fn c Fn is the con- 
served forward axial flux. Conservation of = 
c b Bcb = J2 n I PnC Bn CBn is a trivial issue because cb 
fulfills the linear equation — zcb ~ — Bcb. An analogous 
proof shows that Q F and Qb are also conserved when 
forward components are neglected instead. 

An alternative analysis can be formulated in the light 
of conserved phase symmetries when F — B interactions 
are neglected. Both forward and backward axial fluxes 
can be also envisaged as the conserved charges associ- 
ated to independent global U(l) symmetry transforma- 
tions on F and B components, respectively. This is clear 
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in Ea. 1281 . which is invariant under U(1)f global phase 
transformations 

c F -> e ieF c F (29) 

when M(c F ) = M(e ieF c F ). This is the situation for all 
type of nonlinearities that can be expanded in odd power 
series in conservative lossless systems. Again, the sim- 
plest example is a cubic or Kerr nonlinearity, for which 
M(c) = c*Tc. Since the equation for backward compo- 
nents is basically linear (— its ~ Bc B ), global U(1)b 
invariance (cb — * e j9B ce) is trivially fulfilled as well. 
U(1)f ®U(1)b invariance requires the independent con- 
servation of its associated U(1)f and U(1)b charges, 
which are nothing but the F axial flux (Qf) and the B 
axial flux (—Qb); that is, dQp/dz = and —dQs/dz = 

0. If the situation were the opposite and forward compo- 
nents were neglected, a complete analogous analysis for 
backward components would hold leading to the same 
conclusion. 

When F — B interactions cannot be neglected, the pre- 
vious argument does not hold and neither Qf nor Qb 
are conserved separately in conservative lossless systems. 
From a symmetry point of view, this can be understood 
by the fact that FBEs ijlfijl and ifTTjl are no longer in- 
variant under independent F — B phase transformations 
(c = cp + cb — > c = e l9f cf + e l0B CB) because in that 
case M(c) = M{e i9F c F + c ieB c B ) + M(c F + c B ) = M(c), 
as can be easily checked for a cubic nonlinearity. Us- 
ing the language of group theory, one would state that 
U(1)f <8> U(1)b symmetry is broken. However, there 
is a residual phase symmetry remaining in FBEs iltil 
and i(T7l when U{l)p ® U(1)b invariance is broken by 
F — B interactions. FBEs are still invariant under simul- 
taneous global phase tranformations on F and B com- 
ponents (cf — > c %e cf and cb — > c iS cb, which implies 
c -> e l8 c) provided that M(e l8 c) = M(c). Notice that 
this U(l) symmetry is a particular case of the higher- 
order U(1)f®U(1)b symmetry when Op = and 9b = 0. 
Since U(1)f ® U(1)b symmetry is broken, the U(1)f 
charge, i.e, the F axial flux (Qf), and the U(1)b charge, 

1. e., the B axial flux (—Qb), are not conserved. However, 
in such a situation the remaining U(l) symmetry guaran- 
tees the conservation of a new quantity, namely, the sum 
of the U(1)f and U(1)b charges. That is, Qf + (—Qb) 
has to be conserved. The total axial electromagnetic flux 
Q = Qf — Qb appears then as the conserved charge of 
the [7(1) symmetry associated to the breaking pattern 
U(1)f®U(\) b ^U(1)f+b- 

There are two different and interesting physical sit- 
uations in which the analysis described above for sys- 
tems in which F — B interactions can be neglected is 
valid. The first one is general forward pulse propaga- 
tion, that naturally takes place in the frequency domain 
or, equivalently, in the time domain. The second one is 
monochromatic (or quasi-monochromatic) non-paraxial 
forward evolution describing the nonlinear propagation 
of spatial structures. It is remarkable that despite the 



distinct nature of both phenomena, they are just partic- 
ular cases of FBEs and, thus, equally described by the 
same formalism. We will study them separately in the 
next two subsections. 



A. Forward pulse propagation 

When one neglects F — B interactions, for example 
by eliminating backward components, the most general 
form of FBEs is given by Ea. l28|l . that in continuum 
frequency-notation reads: 

d 

- i—c F (u, z) = f3 n (cj)c F {uJ, z) 
oz 

+ ^rr^Y, I du'M^iu^'-^pyiiu 1 ». (30) 

71' 

The nature of the nonlinearities define the particular 
form of the nonlinear modal matrix function M nn >(w, u'). 
According to what was explained in Section ITTTl the func- 
tional form of the nonlinear function F and, thus, of the 
M matrix function, can be systematically constructed 
out of the mode amplitudes of the linear propagation 
problem together with the standard nonlinear coefficients 
(X (3) ,x (5) , ■•■)• Thus, many formal properties of M will 
be inherited from linear modes. The dependence of M on 
modal indices and frequency will have much to do with 
the particular dependence of the linear mode amplitudes 
on spatial coordinates and frequency. Different physi- 
cally well-grounded assumptions on the properties of M 
on modal (and polarization) indices and frequency can 
be then made by analyzing linear mode characteristics. 
The second element to take into account is the exten- 
sion of the frequency spectrum cf(uj). A considerable 
simplification is achieved for sufficiently narrow spectra, 
whereas wide bandwiths, as those naturally appearing in 
highly-nonlinear fibers (e.g., in supercontinuum genera- 
tion), would demand to consider the frequency depen- 
dence of En. <l30|l in its total extent. 

It is clear that, in the most general case, Eq. jjjDt in- 
volves an intrincate dynamics since the nonlinear matrix 
function M nn '(u),u)') leads to nonzero couplings between 
different frequency, modal and even polarization com- 
ponents (recall that polarization indices are included). 
Thus, even starting from a simple spectral, single-mode, 
single-polarization configuration, if the system had no 
physical mechanism to minimize the great variety of cou- 
plings induced by M, spectral evolution, as described by 
Ea. ^Ojl . would generate a more and more complicated 
spectral function by exciting new frequency, modal and 
polarization components. 

An important feature of Ea. ll3T))l is the existence of 
an inherent spatial modal interference of nonlinear ori- 
gin, represented by the nondiagonal nature of the non- 
linear matrix function in the modal indices -M nn > ^ 
(n ^ n')~ in the most general case. In some special sit- 
uations, however, this modal interplay can be zero or it 
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can be just neglected (M m > = or M nn > « 0, when 
n 7^ n'). Exact cancellation of matrix elements occurs 
by symmetry considerations in most cases. A simple ex- 
ample would be given by a rotationally invariant system. 
Its linear modes are solutions with well-defined angular 
momentum and, consequently, the modal index of the 
spectral function c F is thus labelled by the angular mo- 
mentum index I. In such a case, angular momentum 
conservation requires that no mixing of spectral com- 
ponents with different angular momentum occurs, thus 
eliminating nondiagonal terms in the nonlinear matrix 
function corresponding to different values of the angu- 
lar momentum index I; that is, M has to be diagonal in 
these indices: M nn > ~ 6w. In other cases, some nondi- 
agonal matrix elements can be negligible because of the 
different shapes of linear modes amplitudes involved in 
the calculation of the overlapping integrals appearing in 
the definition of M (for example, for a Kerr cubic nonlin- 
earity: M nn , ~ c* m ,c m J R2 4>* m ,4>* n ,<i) m (f) n ). In some cases, 
these integrals can be very small for reasons that do not 
rely on the presence of particular symmetries. 

Whatever these reasons may be, in the case that modal 
interplay does not exist, or this can be neglected -that 
is, when M nn > = or M nn , » if n ^ n'~, then Eq.ljSOJ 
decouples into independent equations for every mode in- 
dex, 

d 

- i— Cp((J,z) = /3 n (uj)c F {u,z) 
oz 

+ ^TT^ I ^>,w';cf)^(wV)> (31) 

which, in turns, leads to independent conservation laws 
for the different modal F and B axial fluxes: dQp /dz = 
and dQ^/dz = (Q F l) = J^c^cj? 5 and = 

In nonlinear propagation in optical single-mode fibers 
it is commonly assumed that the spatial dependence of 
the electric propagating field is just given by the ampli- 
tude of the fundamental mode of the fiber. In our con- 
text, this statement is equivalent to say that the matrix 
M is one-dimensional and involves the fundamental mode 
only. The equation describing forward nonlinear propa- 
gation in such a fiber would be Ea. ll3"T)) for n = (funda- 
mental mode) . Even if the fiber is not single- mode but in- 
volves modes widely separated by large gaps in /3's, as in 
some highly-nonlinear microestructured fibers, the same 
equation can still remain approximately valid. The rea- 
son is that in such a case intermodal interactions are rel- 
atively suppressed with respect to modal self-interactions 
because of the very different forms of linear mode ampli- 
tudes, originated by both the discrete symmetry of the 
fiber and their very different values of /3, leading to zero 
or small overlapping integrals and, thus, to small val- 
ues of M nn i when n ^ nf. The resulting approximated 
equation — Eq.lplJ restricted to the fundamental mode 
(n = 0) — is equivalent to the Forward Maxwell Equa- 
tion (FME) for a highly-nonlinear microstructured fiber 



used to adequately describe supercontinuum generation 
in this type of fibers Q. 

B. Non-paraxial spatial evolution 

Monochromatic (or, in practice, quasi-monochromatic) 
propagation is also a particular situation that can be de- 
scribed using FBEs. Eqs. ijlfijl and ifTTjl remain cer- 
tainly valid when the frequency content of the F and 
B spectral functions cp and cb is restricted to a sin- 
gle frequency uiq. In such a case, only intermodal inter- 
action or modal self-interaction (including polarization) 
play a role since the nonlinear matrix function has a triv- 
ial dependence on frequency M nn > (u>, u>'; c) — M nn > (loq: c) 
fixed by the propagation frequency ujq. Mathemati- 
cally, FBEs become (in order to symplify the notation: 
Pon = Aj(wo), M nn i(c) = M nn '(u) ;c), c(z) = c(z;lu )) 

d 

- i—c n F {z) = I3 0n c n p{z) 

^POn , 

n' 

d 

-ig-c B (z) = -(3 0n c n B {z) 

^POn , 
n 

and c = cp + cb ■ 

Since the previous FBEs involve spatial modal indices 
only, it is interesting to write them also in the spatial 
domain, so that spatial degrees of freedom appear explic- 
itly. This process is the inverse of the one we followed 
to obtain the second-order modal equation in Section ITT1 
that is, we reintroduce the spatial field amplitudes as 
<f>F,B(x,y,z) = S n c£ )S (z)0 n (z,y), K being the eigen- 
functions of the linear operator Lq = + k^n^ and /3 n o 
their corresponding eigenvalues evaluated at the fixed fre- 
quency ujo — ho/ c. The outcome is 

(~i§- z ~ 4 /2 ) 4>f = \Lv 1/2 M{<i>)4> (34) 

(-^ + 4 /2 ) <t>B = -\L- 1/2 Mm, (35) 

together with cj> = 4>f + 4>b- The nonlinear term M(4>)4> 
includes all type of nonlinearities in the spatial domain 
that permit an expansion in power series. The usual 
example is the Kerr nonlinearity M(<f)<p ~ (4>*4>)4>- 

It is also an interesting exercise to prove that one can 
derive the standard second-order wave equation from the 
first-order spatial FBEs 1-1411 and l|35|) . If we sum and 
substract Eq.<P3J and Eq.JBSJ, we obtain 

-«7tW>f + M-£o /2 (<^-<M = (36) 
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and 

- i^-(fa -fa)- Ll /2 {fa + fa) = Lq 1/2 M(4>)<P, (37) 

respectively. By applying (id/dz) onto Eo. lj36)) . susbti- 
tuting En. lT?7jl into the resulting equation, and taking 
into account that 4> — fa + fa, one gets the Helmholtz- 
type nonlinear wave equation: 

{^+V 2 t +k 2 n 2 +M (0)^ = 0. 

As mentioned in Section IrTfl and it is also evident in 
Ea, lj35ll . the complete elimination of backward compo- 
nents {(f> b — ► 0) implies the vanishing of nonlinearities 
(M — > 0). Conversely, if M ^ then backward compo- 
nents are necessarily generated: fa ^ (even if they did 
not initially exist). Therefore, there exists a close link 
between backward amplitudes and the nonlinear func- 
tion M. Certainly, in most experimental situations in 
which an axially-invariant system is axially illuminated 
along a privileged (say, forward) direction exclusively, 
backward components are small as compared to forward 
amplitudes. However, as backward FBE show, they can- 
not be identically zero in the presence of nonlinearities. 
They would be exactly zero only in the case the sys- 
tem behaved linearly. In such a case, since the system 
is axially invariant and it is illuminated in the forward 
direction only, there would be no axial inhomogeneities 
that could produce linear reflections. Thus, small back- 
ward amplitudes are expected to be generated by small 
nonlinearities so that fa — > when M (linear case) . 
Our interest lies now in the quantification of the rela- 
tion between fa and M in the small backward amplitud 
regime (fa = Sfa < 1). 

In order to clarify the calculation, we parametrize the 
"size" of the nonlinearity by redefining the nonlinear func- 
tion as M — 7M, where 7 is a dimensionless real parame- 
ter. In this way, we can approach the linear case (M — » 0) 
by taking the limit 7 — > 0. It is easy to realize that 7 has 
to be proportional to the input power 7 ~ P when 7 is 
small. Certainly, M is a function of the input power ver- 
ifying that M — > as P — * 0. Assuming, as usual, that 
this dependence is analytical this implies that M ~ P 
when P — > 0. The same argument for 7 leads to M ~ 7 
for small values of 7. Therefore, the small 7 and the 
small P regimes are, in fact, the same one. From a phys- 
ical point of view, this provides a physical meaning to 
the dimensionless parameter 7 in the small nonlinearity 
regime: 7 ~ P. Notice that, in the general case of a 
nonlinearity that can be expanded in power series in <fi, 
the auxiliar nonlinear function M can depend itself on 7. 
However, for our purposes, it is enough to consider that 
M = 0(7) ( and thus M = 0(1)) because we are going 
to be interested in leading order terms. 

Following the reasoning above, we consider the total 
field amplitude as the sum of a large forward amplitude 
and a small backward amplitude: <f> — fa + Sfa. Be- 
sides, the small backward amplitude is a function of the 



nonlinear parameter 7 -Sfa — Sfa(j)- and it has to 
verify that Sfa (0) = 0, since we are assuming that no 
backward radiation is present in the absence of nonlin- 
earities. The nonlinear function can be then expanded as 
M(fa + Sfa) = M(fa) + (8M/d<f>\^) Sfa + O(Sfa) 2 . 
Substituing the total amplitude 4> = fa + Sfa and the 
expansion of M in Eo. ll35)l and keeping the leading order 
terms in Sfa only, one gets 

= ~Lq 1/2 M{fa)fa + 0(Sfa) 2 . 

Since M = 0(7), we can find the leading order term in 
7 for Sfa from the previous equation: 

Sfa = -~ (~i^ z + io /2 ) 1 L Q 1/2 M(fa)fa + 0( 7 ) 2 , 

so that Sfa is also 0(7). We can proceed analogously 
with the forward FBE $M to nnd 

ir^z - L « /2 ) 0f = 

^L- 1/2 M(fa)fa + X -L- 112 (M(fa) + ^F^j Sfa 
+0(Sfa) 2 , 

which, after introducing 7 dependences (M = 0(7), 
Sfa — 0(7)), becomes 

(-i^ - L 1 ^ fa - l -L- 1/2 M{fa)fa + 0( 7 ) 2 . (38) 

Therefore, neglecting backward components is equiv- 
alent to consider FBEs up to 0(^) 2 terms. Or, equiv- 
alent^, the pure forward equation (that is, FBEs with 
4> ~ fa and fa w 0) is just a weak-nonlinearity approx- 
imation. More specificaly, it is the leading-order contri- 
bution in the nonlinear parameter 7 to FBEs. For this 
reason, it is possible to find equivalent forms to the for- 
ward equation different than l(38jl . An interesting alter- 
native version of the forward equation is easily obtained 
by using the following property 

L\ /2 fa + h 1/2 M(fa)fa + 0( 7 ) 2 = 

Lo 2 (l + l^Mifa)) fa + 0( 7 ) 2 = 
(L Q + M( fa)) 1/2 fa + 0(j) 2 , 
which allows us to write Eq.JHEJ as 

- i^-fa = (L + M{fa)) 1/2 fa + 0( 7 ) 2 . (39) 

The previous version of the forward equation has been 
used to simulate monochromatic nonlinear propagation 



10 



of spatial structures in photonic crystal fibers m the 
non-paraxial regime. Despite the pure forward equation 
is first-order in z, it has an intrinsic non-paraxial nature. 
Unlike in the nonlinear Schrodinger equation, the stan- 
dard integral J 4>* f 4>f is not the conserved quantity. This 
can be clearly seen if one writes the evolution operator 
associated to Eo. ll39)) for an infinitesimal axial step e: 

1/2 

(j) F (z + e) = cxpie (L + M(4> F {z))) 4>f(z). 

This evolution operator is not unitary because the op- 
erator Lq + M(4>f(z)), despite it is self-adjoint, is not 
positive-definite, inasmuch as it can have negative eigen- 
values corresponding to evanescent waves (/3 2 < 0). The 
loss of unitarity is due to this evanescent modes leading 
to the non-conservation of the integral J (jf F 4>F- 

VII. CONCLUSIONS 

The experimental availability of high nonlinearities is 
expected to unveil a number of new effects that will 
force us to extract information from Maxwell's equa- 
tions in a more accurate manner. In this paper we shed 
some light into the problem of revealing the close inter- 
play between backward components and nonlinearities 
in axially-invariant systems. With a minimum amount of 
approximations, we have been able to find a system of two 
coupled first-order equations for the forward and back- 
ward spectral components of the electromagnetic field, 
the so-called forward-backward equations. The explicit 
appearence of forward and backward components as well 
as of their nonlinear couplings in these equations is useful 
to quantify under which conditions nonlinearly-generated 
backward components can be relevant in a new scenario 
of highly-nonlinear effects. 

From the formal point of view, the FBEs are spe- 
cially appealing in the sense that they admit a simple 
bi-spinor representation that closely resemble that of a 
Dirac equation for the positive and negative components 
of a fermion wave function in 1+1 dimensions. Similarly 
to Dirac equation, the use of the algebraic properties of 
the spinor FBEs allows us to obtain the conserved quan- 
tities associated to them in an elegant way. In the same 
manner, all conserved quantities also admit an interpre- 
tation as conserved charges associated to phase symme- 
tries. 

The dimensional reduction is a remarkable issue of 
the modal approach followed here. The original 3+1 di- 
mensions (3 spatial, 1 frequency) of the starting wave- 
equation for E w (x, y, z) (Eq.JJJl) are reduced to 1+1 (1 
spatial, 1 frequency) in the FBEs. The modal approach is 
a way of "integrating out" the transverse spatial degrees 
of freedom (x and y). Of course, the coupling between 
tranverse degrees of freedom in Eq.QJ does not dissa- 
pear in FBEs. It transforms into the couplings between 
the different modal components of the spectral function 
c n which are mathematically encoded in the nonlinear 



matrix function M nn i. For the case in which only a 
few modal components are relevant, the process of di- 
mensional reduction provides a dramatical simplification. 
Typical envelope equations for propagating pulses are the 
result of a similar process in which the propagation of 
only one linear mode (usually, the fundamental mode) is 
assumed. The FBEs, however, provide a natural way of 
dealing with a more complex modal structure and, more- 
over, they allow one to directly work with the frequency 
content of the propagating field; that is, with its spectral 
components c n (cu, z). In this sense, there is no need to re- 
sort to the concept of pulse envelope (unless one wants to 
make contact to other approaches). For the same reason, 
the FBEs are equally suitable to describe pulse propaga- 
tion with extremely large bandwiths since no assumption 
on the form of c„(w, z) is required. Ideally, they could 
handle any type of temporal-spectral behavior provided 
the only assumption needed, the "scalar approximation" 
(in its strong or weak form — see section HT1— ). is reason- 
ably fulfilled. 

Another scenario in which these equations can be use- 
ful is that in which there are relevant spatial-temporal 
effects. The FBEs inherently include couplings between 
frequency and modal indices through the nonlinear ma- 
trix function M nn /(w,u/; c). Since modal indices corre- 
spond to spatial amplitudes of linear modes, the nonlin- 
ear matrix function M nn > (u>, u>'; c) is constructed out of 
these amplitudes as an overlapping integral (section ITU . 
Thus, part of the frequency dependence of M nn > is due 
to the explicit dependence of linear mode amplitudes on 
frequency. In the ultra-wide spectrum regime, this de- 
pendence cannot be neglected and, therefore, the contri- 
bution of several spatial modes can produce simultaneous 
couplings between modal indices and frequencies which 
can be naturally treated in the framework of the FBEs. 
In terms of the original 3+1 Maxwell's equations, these 
effects would correspond to spatial-temporal phenomena, 
in which the spatial and temporal degrees of freedom of 
the electric field E(x,y,x,t) could not be factorized. 

Sumarizing, the FBEs provide a general framework 
to deal with nonlinearities in axially-invariant inhomoge- 
nous dielectric media, limited only to a reasonable valid- 
ity of the "scalar approximation" (in its strong or weak 
form) . As we have seen in section lVll its generality can be 
made evident after analizing two apparently disconnected 
cases: forward pulse propagation (a purely temporal phe- 
nomenon) and monochromatic non-paraxial evolution of 
spatial structures (a purely spatial phenomenon). Both 
are equally and naturally described within the FBEs for- 
malism. The fast progress in nonlinear optics experi- 
ments provides effects of increasing complexity both in 
the spatial and time domains as well as in the interplay 
between forward and backward components. Strong cor- 
relations between spatial and time domains and forward 
and backward components will play a more and more im- 
portant role. In this context, the FBEs can be a suitable 
and convenient tool to encompass a variety of different 
phenomena within a common framework. 
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